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Quantum Monte Carlo calculation of Compton profiles of solid lithium 
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Recent high resolution Compton scattering experiments in lithium have shown significant discrep- 
ancies with conventional band theoretical results. We present a pseudopotential quantum Monte 
Carlo study of electron-electron and electron-ion correlation effects on the momentum distribu- 
tion of lithium. We compute the correlation correction to the valence Compton profiles obtained 
within Kohn-Sham density functional theory in the local density approximation and determine that 
electronic correlation does not account for the discrepancy with the experimental results. Our calcu- 
lations lead do different conclusions than recent GW studies and indicate that other effects (thermal 
disorder, core- valence separation etc.) must be invoked to explain the discrepancy with experiments. 



I. INTRODUCTION 

Inelastic x-ray scattering is called Compton scattering 
when the energy and the momentum transferred are large 
compared to the characteristic energy and reciprocal in- 
terelectronic distance of the scattering system. Comp- 
ton scattering probes the electronic structure of materi- 
als through the electronic momentum distribution and, if 
the scattering system is metallic, gives direct information 
on various characteristics of the Fermi surface, such as 
position and size of the Fermi breaks and their renormal- 
ization due to electron-electron correlation. Fermi mo- 
menta are of the order of 1 a.u. so that high-resolution 
in momentum are necessary in order to resolve features 
related to the Fermi surface. In the last few years, the 
advent of high intensity, high energy and well-polarized 
synchrotron sources has made possible to obtain resolu- 
tions of the order of 0.1 a.u. and high statistics. Record 
resolutions of the order of 0.02 a.u. have been recently 
achieved for solid beryllium 0] and lithium 0]. 

In the range of energy and momentum transferred where 
the recoil electron can be considered free, i.e. within 
the impulse approximation ||, the experimental double- 
differential Compton cross section is related to the mo- 
mentum distribution n(p) of the electronic system as (in 
atomic units K = e = m = 1): 
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(1) 



where lu and q are the energy and momentum transferred, 
loi is the energy of the scattered photon and (dcr/df2) 
is the Thomson differential cross section. The outcome 
of a Compton scattering experiment is therefore given by 
what is called a Compton profile in a given direction q 
and defined as 



J(p) = / dpn(p)<J(p-q-p). 



(2) 



The Compton profile at p is the integral of the momen- 
tum distribution on a plane perpendicular to the unit 
vector q at a distance p from the origin. The momentum 
distribution is expressed in terms of the wave function of 
the electronic system as 

"(P) = y[dri...dr N J dr' e 4pr ' **(r x , . . . ,r N ) x 

9( ri +i>,..., rff ), (3) 



where N is the number of electrons and V the volume of 
the system. 

Since Compton experiments can only access the momen- 
tum distribution in an indirect way and there are diffi- 
culties in handling background or subtracting the core 
contribution (the impulse approximation may not hold 
for the core electrons it is important to estab- 

lish the performance of the technique for a system where 
high resolution experimental data are available and cor- 
related calculations can be performed. Being a low-Z 
material, lithium has been the subject of Compton stud- 
ies since the early days of x-ray Compton scattering || . 
Recent high-resolution experiments have been conducted 
on lithium and Fermi surface signatures have been inves- 
tigated either by analyzing first and second derivatives of 
the Compton profiles JtJ or attempting a reconstruction 
of the momentum distribution Q . 

Despite overall shape similarities, there is a clear discrep- 
ancy between the measured Compton profiles of lithium 
and the theoretical profiles computed within Kohn-Sham 
density functional theory in the local density approx- 
imation (LDA) even when the Lam-Platzman correla- 
tion corrections M are included. These discrepancies 
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were attributed to inadequate treatment of correlation 
in the Kohn-Sham single particle picture and it was 
proposed that plasmaron (a bound hole-plasmon state) 
losses could explain the observed difference ||. Two re- 
cent theoretical studies have tried to explain this dis- 
crepancy but reached opposite conclusions. Kubo calcu- 
lated the Compton profiles using the occupation num- 
bers obtained from a GW calculation and recovered a 
good agreement with experimental results p0| . Dung- 
dale and Jarlborg simulated the effect of thermal disorder 
on the Kohn-Sham Compton profiles, found that disor- 
der leads to a derealization of the momentum distribu- 
tion and concluded that this effect, combined with the 
Lam-Platzman correlation corrections, accounts for the 
discrepancy with the experimental Compton profiles [[llj . 

In order to determine the correction to the LDA Compton 
profiles due to electron-electron correlation, we perform 
a fully correlated calculation of the momentum distri- 
bution of solid lithium within pseudopotential quantum 
Monte Carlo (QMC). The QMC Compton profiles differ 
from the experimental data by Sakurai et al. and in- 
dicate that Lam-Platzman corrections give a satisfactory 
description of electronic correlation in this system even 
though they are isotropic and cannot reproduce the ob- 
served directional dependence of the QMC corrections to 
the LDA Compton profiles. Therefore, the QMC results 
for bec lithium differ from the GW calculations |l(J and 
suggest that the discrepancy between conventional band 
theory and experiments originates from other sources: 
temperature effects JTTJ or problems in the interpretation 
of the experimental data are possible explanations. 

Recently, QMC calculations of Compton profiles have 
also been performed for a rather different system, sili- 
con jl^]. Due to the absence of a Fermi surface, the mo- 
mentum distribution of silicon is a smooth function and 
finite size errors in QMC are easier to correct than in a 
metal. These studies concluded that, also in silicon, cor- 
relations effects are well described by the Lam-Platzman 
corrections and do not fully account for the discrepancy 
between the theroretical Kohn-Sham LDA profiles and 
the experimental results. 

In Section |l[ the characteristics of our density functional 
theory calculations are outlined. In Section III, the func- 
tional form of the QMC wave function is described. In 



Section [V, we compare the Kohn-Sham Compton pro- 



files obtained from pseudopotential and full-core LDA 
calculations. We then discuss the role of electron- ion and 
electron-electron correlation on the momentum distribu- 
tion in QMC. We finally present the QMC correlation 
corrections to the valence Kohn-Sham Compton profiles 
and compare them with the experimental results. In Ap- 
pendix |A|, the linear tetrahedron method to obtain the 
Kohn-Sham Compton profiles is briefly outlined. In Ap- 
pendix [§], convergence of valence properties of lithium 



from full-core plane-wave calculations is discussed. Vari- 
ance minimization, variational Monte Carlo (VMC) and 
diffusion Monte Carlo (DMC) methods are briefly pre- 
sented in Appendix 



II. LDA CALCULATIONS 

We study bec lithium at the experimental lattice constant 
of 6.60 a.u. We carry out pseudopotential and full-core 
calculations within LDA density functional theory in a 
plane-wave basis. The pseudopotential calculations pro- 
vide the single-particle orbitals that enter in the QMC 
wave function and the reference momentum distribution 
for the correlation corrections determined within QMC. 
The full-core calculations of the valence contribution to 
the Compton profiles are carried out to account prop- 
erly for core- valence orthogonality. In Appendix |X[ we 
describe how to evaluate the Kohn-Sham momentum dis- 
tribution and construct the Compton profiles using the 
linear tetrahedron method. 

In the pseudopotential calculations, we use the Troullier- 
Martins pseudopotentials whose s and p components are 
generated with a cutoff radius of 2.4 a.u. Jl3|]. The plane- 
wave cutoff is set to 16 Ry and 44 special fc-points H] in 
the irreducible wedge of the zone are used for zone sam- 
pling during iteration to self-consistency. The Compton 
profiles are generated from 16206 fc-points in the irre- 
ducible zone unfolded in a sphere of 2 a.u. radius (the 
mesh spacing is 0.0136 a.u.). The profiles arc converged 
with respect to mesh spacing and sphere radius but ap- 
pear to be very sensitive to the value of the Fermi energy 
that must be carefully determined as described in Ap- 
pendix |A|. 

The valence Compton profiles from a plane- wave full-core 
calculation are obtained with a plane-wave cutoff of 400 
Ry and 3311 fc-points in the irreducible zone, unfolded 
into a sphere of 4 a.u. radius (the mesh spacing is 0.0238 
a.u.). The Compton profiles arc in good agreement with 
the profiles computed with the linear augmented plane 
wave method [ fl5| |. In Appendix [b|, we discuss conver- 
gence issues for valence properties computed from a full- 
core calculation in a plane- wave basis. 



III. QUANTUM MONTE CARLO 
CALCULATIONS 



In the quantum Monte Carlo simulations of bec lithium, 
we treat lithium as a pseudo-ion of charge one plus one 
valence electron. We test the use of both a local and a 
non-local pseudopotential Jl6[ |. To simulate an infinite 
solid, we model the system as a collection of ions and 
electrons in a simulation cell periodically repeated. An 
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effective computation of the Ewald sums over the im- 
ages of the potential is obtained by optimizing the sep- 
aration between the long- and short-range components 
of the electron-electron, electron-ion and ion-ion interac- 
tions We will consider cubic simulations cells con- 
taining 54, 250 and 686 lithium atoms. 

The wave function used in these calculations is a deter- 
minant of single-particle orbitals multiplied by a Jastrow 
factor describing electron-electron and electron-ion cor- 
relations: 



$ = D f xD 1 exp 



JV 



JV 



(4) 



stars, x s , are optimized using the variance minimization 
method fHH. 

A brief description of variance minimization, VMC and 
DMC methods is given in Appendix ^|. 



A. fe-point sampling 

We impose that the wave function satisfies boundary con- 
ditions that can be either periodic or arbitrary: 



• ,rj 



,rjv) = 
e ik » R »*(ri 



-,t n ), (9) 



DT and D x are the Slater determinants of single parti- 
cle orbitals for the up and down electrons respectively, u 
correlates pairs of electrons and \ is a single-body term. 

We consider two forms of single-particle orbitals. 1) We 
test the very simple choice: 



0o (r) 



(5) 



2) We compute the single-particle orbitals from a LDA 
density functional theory calculation within a plane-wave 
basis: 



'LDA 



£ Ck+ G^ k+G >", 

G 



(6) 



where G are the reciprocal lattice vector of the underly- 
ing bcc lattice. We discuss the occupation of the orbitals 
below. 

The electron-electron term in the Jastrow factor is peri- 
odic over the cell and contains no free parameters. It is 
derived within the random phase approximation (RPA) 
and has been extensively used for the homogeneous elec- 
tron gas [jOI and solid hydrogen calculations |filj|| . It de- 
scribes both the exact short range (cusp condition) and 
large distance (plasmon) behavior. We only impose the 
antiparallel cusp conditions: 



Mi 



(ry) = M TT (ry) = u ll (r lj ) = urpa^). 



(7) 



The single-body term in the Jastrow factor has the peri- 
odicity of the underlying lattice, is expanded in Fourier 
components and rewritten as the sum over stars of recip- 
rocal lattice vectors: 



(8) 



GEs 



Since the point group of bcc lithium is symmorphic, 
the phases $g can be set to unity. We include up to 
twelve stars but convergence is obtained already with 
an expansion over seven stars. The coefficients of the 



where one electron has been displaced by a translational 
vector of the simulation cell, R s . The Jastrow component 
is periodic over the cell and does not affect the phase of 
the wave function. The determinant satisfies the above 
equation if the single-particle orbitals, <p, obey similar 
equations, 4>(r + R s ) = exp{ik s • R s }0(r), that is if 
the orbitals correspond to a single fc-point sampling of 
the simulation cell given by k s . If we impose periodic 
boundary conditions on our cubic cell of side L, the or- 
bitals are Bloch states with wave vectors defined on a 
cubic grid centered on the origin with spacing 2ir/L. For 
arbitrary boundary conditions, the grid of wave vectors 
is shifted by k s . 

The idea of /c-point sampling in QMC was introduced to 
achieve a faster convergence in the total energy versus 
system size, possibly faster than in a periodic calcula- 
tion p2| |. For bcc lithium, since the valence band is an 
s-band, a calculation not at the T-point yields a higher 
energy and corresponds to an excited state of our finite 
simulation cell. On the other hand, the allowed momenta 
coincide with the wave vectors compatible with the condi- 
tion of periodicity on the simulation cell, so we can obtain 
a higher resolution in momentum space by performing 
calculations with different boundary conditions Jl2| . 

We only consider /c-samplings that yield a grid of fc-points 
with inversion symmetry so that we can construct a real 
wave function by occupying linear combinations of pairs 
of orbitals: 



1 



1 

2i { 



c). 



(10) 



By restricting ourselves to real wave functions, we can 
perform both VMC and DMC calculations and avoid the 
complications of dealing with complex wave functions in 
DMC J23]. For a cubic simulation cell, there are only 
four vectors k s that preserve inversion symmetry corre- 
sponding to r, X, M and R sampling of the cubic simu- 
lation cell. The T-point calculation yields a cubic grid of 
k- vectors centered on the origin with spacing 2tt/L and 
the additional calculations provide three grids shifted by 
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(2tt/L)/2 in the [100], [110] and [111] directions respec- 
tively. 

To construct the determinantal part of the wave func- 
tion, we compute the orbitals on the grid of fc-vectors 
compatible with the boundary conditions and occupy the 
orbitals, 4>q or 0lda> within the first Brillouin zone with 
the lowest free-electron (<fro) or LDA (0lda) energy. In 
general, the highest occupied level is degenerate and only 
partially occupied, so we should employ a linear combi- 
nation of determinants in order to construct a wave func- 
tion with the proper symmetry. We instead always use a 
single determinant and symmetrize the momentum dis- 
tribution on the grid of allowed wave vectors by averaging 
it over symmetry related A:- vectors on the grid. This pro- 
cedure is done separately for the four fc-point samplings. 
Finally, by combining the results of the four calculations 
and applying the symmetry rotations of the cubic group, 
we obtain a momentum distribution defined on a mesh 
with spacing (2tt/L)/2. 

We consider cubic cells with 54, 250 and 686 atoms. For 
the 54 atom cell, we only carry out a T-point calculation. 
The wave function is constructed from the orbitals corre- 
sponding to the lowest four complete shells of fc-vectors, 
so it is real, has the full symmetry of the lattice and is a 
spin-singlet. For the 250 and 686 atom cell, we compute 
the momentum distribution for four wave functions cor- 
responding to r, X, M and R sampling of the cubic cell. 
The grid in momentum space has spacing 0.095 a.u. and 
0.068 a.u. for the 250 and the 686 atom cell respectively. 



B. Momentum distribution and Compton profiles 

The momentum distribution of the correlated wave func- 
tion is computed as the expectation value over the distri- 
bution given either by 'J 2 (VMC) or by the product of the 
trial wave function and the fixed-node solution (DMC): 



n(p) = 

y- 



dr 



, riv . r ' i H{r 1 ,..,ri-\-r',..,v N ) 



(11) 



At a given Monte Carlo step, we uniformly sample M 
random positions, r', within the simulation cell and, for 
each position and particle, compute in turn the above 
ratio [ p4| . The value of M depends on the size of the 
system and is determined to optimize the efficiency in 
sampling the momentum distribution. The allowed val- 
ues of momentum coincide with the sets of vectors com- 
patible with the condition of periodicity imposed on the 
simulation cell (Eq. ||). 

To obtain the Compton profiles, we have to integrate the 
momentum distribution over planes perpendicular to a 
given direction. It is possible to directly compute the 



Compton profiles within QMC. For instance, to evaluate 
the profile in the [100] direction, we could estimate the 
following expectation value: 



jfe) = E 



(27T) 
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(12) 



This procedure is equivalent to evaluating the integral 
of the momentum distribution as the histogram over the 
momenta compatible with the boundary conditions: 



J(p x )^(2ir/L) 2 n (P)' 



(13) 



and it is clearly a poor representation of the integral of a 
function with several discontinuities, especially if the grid 
of fc-vectors is coarse. The Compton profiles computed 
in this way show indeed a strong dependence on the size 
of the simulation cell. 

To obtain Compton profiles with reduced finite size er- 
rors, we want to integrate a smoother function than 
the QMC momentum distribution, jiqmc(p)- We have 
to select a reference distribution, «model(p), whose 
Compton profiles, Jmodel(p), can be computed with 
high accuracy and such that the difference, An(p) = 
n QMc(p) — "-model (p) 7 is a smoother function than the 
original Uqmc(p)- The Compton profiles are then ob- 
tained as 



Jqmc(p) = ./model (p) + AJ(p) 



(14) 



where the corrections, AJ(p), is computed by integrating 
An(p) on the grid defined by the four fc-point samplings 
using the linear tetrahedron method. The only difference 
in computing AJ(p) from the method described in Ap- 
pendix |A| is that there is no Kohn-Sham energy defined 
on the grid and no Fermi energy but all grid points are 
considered for integration. In Section IV, we will show 
our choice for the reference momentum distribution. 



IV. RESULTS AND DISCUSSION 



We present calculations of the momentum distribution 
and Compton profiles within LDA Kohn-Sham density 
functional theory. We perform both pseudopotential and 
full-core calculations. We then compute the momentum 
distribution within pseudopotential QMC and determine 
the correlation contribution to the directional Compton 
profiles. 
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A. Kohn-Sham LDA Compton profiles 

In Fig. |], we show the LDA momentum distribution along 
the [110] and [111] directions for a local (s-component) 
and a non-local (s- and p-components) pseudopotcntial. 
The effect of the lattice on the momentum distribution is 
more pronounced when using a non-local pseudopoten- 
tial, yielding a bigger reduction below the Fermi break 
and a more significant contribution from the secondary 
Fermi surfaces (unklapp processes). For both the local 
and non-local pseudopotentials, the momentum distribu- 
tion is strongly anisotropic and goes to zero after the 
first Fermi break since the Fermi surface of lithium is 
completely contained within the first Brillouin zone. 
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FIG. 1. LDA momentum distribution in the [110] and [111] 
directions using a local and a non-local pseudopotential. 

We find that the momentum distribution obtained with 
the non-local pseudopotential is in better qualitative 
agreement with our full-core calculations. Therefore, for 
a more realistic description of bec lithium, we have to use 
a non-local pseudopotential even though this is compu- 
tationally more demanding in QMC and introduces the 
additional locality approximation in DMC 

The agreement of the valence momentum distributions 
from a non-local pseudopotential and a full-core calcu- 
lation is only qualitative because of the pseudopotential 
approximation (lack of correct orthogonalization between 
core and valence). When computing the momentum dis- 
tribution within a pseudopotential scheme, we are under- 
estimating the momentum distribution at high momenta 
and, consequently, overestimating it at low momenta. 
This effect is also evident in the Compton profiles. 

In Fig ^|, the pseudopotential Compton profiles are com- 
pared with the valence profiles obtained with the full-core 



potential in the [100], [110] and [111] directions. We nor- 
malize the Compton profiles to the number of valence 
electrons per unit cell, in this case one electron. As ex- 
pected, the Compton profiles constructed from valence 
orbitals correctly orthogonalized to the core orbitals show 
a higher tail at high momenta than the pseudopotential 
profiles and, since they must integrate to the same value, 
a significantly lower value at low momenta. The Comp- 
ton profile of a free particle system at the same density 
of lithium (r s — 3.25 a.u.) is an up-side-down parabola 
terminating at p F = 0.5905 a.u. Due to electron-ion cor- 
relation, the Kohn-Sham profiles develop a tail beyond 
p F but the sharp features of the presence of a Fermi sur- 
face is clearly visible in any direction at about 0.6 a.u. Its 
location varies however for the different directions as a re- 
sult of the already observed anisotropy (see Fig. |l|) . The 
effect of secondary Fermi surfaces (unklapp processes) 
appears both in the pseudopotential and full-core poten- 
tial valence profiles. The effect is more pronounced in 
the [110] direction at about 0.8 a.u. and all profiles show 
an additional bump in the far tail (above 1.6 a.u.). 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 
p (a.u.) 



FIG. 2. Valence Compton profiles of Li in the [100], [110] 
and [111] directions. The LDA Compton profiles constructed 
using the Troullier-Martins pseudopotential are compared 
with the LDA profiles obtained with the full-core — 3/r po- 
tential. 

Within pseudopotential quantum Monte Carlo, we will 
compute a correction to the Compton profiles due to elec- 
tronic correlations. This correction will be summed to 
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the valence Compton profile from a full-core calculation 
to account properly for core- valence orthogonality. 



B. Quantum Monte Carlo Compton profiles 

In order to calculate the correlation corrections to the 
Kohn-Sham Compton profiles of Fig. |2|, we first deter- 
mine potential and wave function needed to correctly de- 
scribe the momentum distribution of solid lithium within 
QMC. We start modeling bcc lithium with the simplest 
potential and wave function and improve upon it with a 
more sophisticated wave function and Hamiltonian. 

The tests are conducted on a small simulation cell of 54 
electrons and 54 lithium ions on a bcc lattice and with 
periodic boundary conditions. We carry out QMC calcu- 
lations using the Troullier-Martins pseudopotential and 
test the use of a local (s-component) and a non-local 
(s- and p-components) pseudopotential as we have done 
within LDA density functional theory. For either poten- 
tial, we employ both the free-electron orbitals 0o (Eq. ||) 
and LDA orbitals (Eq. ^) determined with a plane-wave 
cutoff of 16 Ry. The Jastrow factor is separately opti- 
mized in each case using 2000 configurations within vari- 
ance minimization. The VMC and DMC energies ob- 
tained with the non-local pseudopotential and the LDA 
orbitals are given by Evmc = —0.2524(1) Hartree and 
Edmc = -0.2591(1) Hartree. 

In Fig |^, we show the spherical average of the momen- 
tum distribution in the case of local and non-local pseu- 
dopotentials. The VMC and DMC spherical momentum 
distribution for bcc lithium are compared with the VMC 
momentum distribution of the homogeneous electron gas 
at the same density (r s = 3.25 a.u.). The wave function 
for the electron gas is given by the product of a deter- 
minant of simple plane waves, <po, and a Jastrow factor 
only containing the electron-electron term described in 
Section III. The momentum distribution of the electron 
gas shows a discontinuity at p F . It is reduced below p F 
and develops a tail at high momenta with respect to the 
non-interacting step- function. The persistence af a Fermi 
break and its location are consistent with the Fermi liquid 
behavior of the system. In the presence of an electron- 
ion pseudopotential, we have a further reduction of the 
size of the Fermi break due to electron-ion correlation. 
In all cases, the VMC and DMC momentum distribution 
are almost indistinguishable indicating that for this sys- 
tem the variational wave function is quite close to the 
fixed-node solution in describing this property. 

For both local and non-local pseudopotentials, we test 
the use of free-electron orbitals, 4>o, versus LDA orbitals. 
The use of orbitals given by a single plane-wave is jus- 
tified if the dependence of the periodic component of 
the Bloch states is not strongly varying with k. This 



dependence can then be simply included in the Jastrow 
part as an electron-ion term, \ (Eq- |^)- As shown in 
Fig. ||, this argument holds in the case of a local po- 
tential: the momentum distribution is not significantly 
different when using the orbitals <fio or the LDA orbitals. 
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FIG. 3. VMC and DMC spherical momentum distribution 
for bcc lithium with a local (upper plot) and non-local (lower 
plot) pseudopotential and a simulation cell of 54 atoms with 
periodic boundary conditions. Either simple plane-wave (4>o) 
or LDA orbitals (<^lda) are employed in the determinantal 
part of the wave function. The VMC momentum distribution 
for the electron gas is also shown. The statistical errors are 
smaller than the size of the symbols. The vertical dashed line 
indicates the position of pf- 
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On the other hand, free-electron orbitals in the determi- 
nant plus an electron-ion term in the Jastrow component 
give a poor description of electron-ion correlation when 
using a non-local pseudopotential. LDA orbitals yield 
a larger reduction below the Fermi wave vector and the 
appearance of additional structure in the tail of the 
momentum distribution as contribution from secondary 
Fermi surfaces. This effect is smeared out because of 
spherical averaging and is more evident when plotting 
the momentum distribution along different directions. 

Having identified the necessary features to describe bcc 
lithium (a non-local pseudopotential and LDA orbitals in 
the determinatal part of the wave function), we perform 
QMC calculations for larger systems and different bound- 
ary conditions to obtain a higher resolution in momentum 
space. We consider two systems with 250 and 686 atoms 
and compute the VMC momentum distribution for four 
wave functions corresponding to T, X, M and R sampling 
of the cubic cell. For the system with 250 atoms, we also 
determine the momentum distribution within DMC. The 
parameters used in the electron-ion Jastrow component 
of the four wave functions for both system sizes are the 
ones optimized for the 54 atom simulation cell. 

As already explained, each fc-point sampling defines a 
grid of allowed wave vectors and the QMC wave function 
is constructed from the orbitals with the lowest LDA en- 
ergy. This procedure yields a different LDA Fermi energy, 
£p a , and set of occupation numbers, for each /c-point 
sampling. The occupation /£ s is equal to two below ep s 
and in general fractional at the Fermi level. In the fol- 
lowing, when comparing the LDA and QMC momentum 
distribution on the grid, we are referring to 



^lda(p) = 



(e(p) 



4 s ), 



(15) 



where c p is the Fourier component of an LDA orbital 
and the momentum p is on the grid corresponding to 
k s -sampling. This ensures that both the LDA and 
the QMC momentum distribution satisfy the sum sule 
J2i n (Pi) = N for each fc-point sampling. 

In Fig. ||, we plot the VMC and DMC momentum distri- 
butions in the [100], [110] and [111] directions for the 250 
atom cell. The LDA distribution evaluated on the same 
fc-vectors is also shown for comparison. The VMC and 
DMC results are quite close to each other with the DMC 
momentum distribution being slightly higher at low mo- 
menta and lower at high momenta than the VMC one. In 
the QMC distribution, we observe a reduction of the LDA 
momentum distribution below the Fermi wave vector, the 
persistence of the discontinuity and an enhancement at 
high momenta due to electron-electron correlation. The 
momentum distribution for lithium both in LDA and in 
QMC is strongly anisotropic. 
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FIG. 4. VMC and DMC momentum distribution of bcc 
lithium in the [100], [110] and [111] directions. The statistical 
error is smaller than the size of the symbols. The momen- 
tum distribution is computed for a 250 atom cell and four 
different fc-point samplings. The LDA momentum distribu- 
tion evaluated at the same vectors is shown for comparison. 



As explained in Section III, we want to define a reference 
momentum distribution such that its difference with the 
QMC distribution is a smoother function to integrate and 
Compton profiles with reduced finite size errors can be 
obtained. In Fig. [|, we show the difference of the QMC 
and LDA momentum distributions along the [110] direc- 
tion for the 250 atom cell. This difference has smaller 
discontinuity than the original VMC or DMC data and 
its integral is less sensitive to finite size errors. However, 
a yet better choice is to compute the following difference: 



An Q (p) = tiqmc(p) ~ ev Tilda (p), 



(16) 



where the parameter a is choosen to reduce the size of 
the discontinuity in the function An". To determine a, 
we minimize the cost function: 

3 



E t A n " (p) - A n ° (p + d ^ )] 2 "qmc (p) 5 



(17) 



where dpi is one grid spacing in the x, y and z directions. 

For both the 250 and the 686 atom simulation cell, we 
find the optimal value of a = 0.8. In Fig. ||, An a (p) 
is plotted in the [110] direction for the 250 atom sys- 
tem. An"(p) is a more regular function to integrate 
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with no significant discontinuity and a similar behavior 
is observed in the [100] and [111] directions. We show be- 
low that this procedure of defining a reference momentum 
distribution is successful since it yields Compton profiles 
that lie on the same curve for the various system sizes 
(see Figs. || and ^). 
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FIG. 5. Difference Bvmc — w LD a and n V MC — a^LDA of 
the VMC and LDA momentum distribution. The momen- 
tum distribution is computed for a 250 atom cell and four 
different fc-point samplings. The parameter a is equal to 0.8. 

To obtain the Compton profiles, we compute A J a (p) by 
integrating An a (p) with the linear tetrahedron method. 
The momentum space is divided in cubes, the cubes in 
tctrahcdra and a contribution to the integral from each 
tetrahedron is computed as in Appendix |a|. A J a (p) rep- 
resents the correlation contribution to the Compton pro- 
file that is finally obtained as 

Jqmc(p) = a JldaOp) + A J a (p), (18) 

where Jlda(p) is the LDA Compton profile. To prop- 
erly take in account core-valence orthogonality, we sum 
the correlation contribution A J a (p) to the LDA valence 
Compton profile obtained from a full-core calculation (see 
Fig. ||) . A J™ (p) is computed using a pseudopotential 
but, since it is the difference of results from two pseu- 
dopotential calculations (LDA and QMC), there is some 
degree of cancellation in the pseudopotential error. 

In Fig. |, the VMC Compton profiles in the [100], [110] 
and [111] directions for the 250 and 686 atom cell are 
compared with the LDA valence profiles from the full- 
core calculation and the experimental results by Saku- 
rai et al. M. The agreement between the simulations 
with 250 and 686 electrons implies that the differences 
between the QMC profiles and either the LDA or the ex- 
perimental curves are beyond finite size errors. Due to 
electron-electron interaction, the QMC profiles are lower 



at low momenta and higher at intermediate and high mo- 
menta than the LDA profiles. The experimental data on 
the other hand show opposite trends. They are lower 
than the QMC profiles below the Fermi wave vector and 
significantly higher in the intermediate and far tail. We 
do not convolute our results with the experimental reso- 
lution since the momentum resolution for the 250 atom 
simulation is only marginally better than the experimen- 
tal resolution of .12 a.u. The continuum line for the QMC 
profiles is a result of the linear tetrahedron construction 
on An Q (p). 

1.3 i 1 1 1 1 1 1 1 1 1 1 1.3 



QMC valence Compton profiles - 1.2 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 
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FIG. 6. Valence Compton profiles for lithium in the [100], 
[110] and [111] directions. The VMC results for the cell with 
250 and 686 atoms and four fc-point sampling are compared 
with the LDA valence profiles from a full-core calculation and 
the experimental results [Q|. 

In Fig0, we show the difference of the QMC and LDA va- 
lence Compton profiles, Jqmc — Jlda, in the [100], [110] 
and [111] directions, so details in the correlation correc- 
tions can be better appreciated. We plot the VMC and 
DMC results for the 250 atom cell and the VMC results 
for the 686 atom cell. The agreement between the results 
of the simulations with the two system sizes is quite good 
with the exception of the [100] direction at low momenta 
as a result of finite size errors for the 250 atom cell. As al- 
ready observed for the momentum distribution, the DMC 
results are only slightly higher at low momenta and lower 
at high momenta than the VMC ones. The difference of 
the LDA and QMC profiles is anisotropic. In the [100] 
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and [111] directions, the QMC profiles do not show the 
sharp features of the LDA profiles at the Fermi break but 
are more rounded off and, consequently, their difference 
is peaked at the Fermi break. 
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FIG. 7. Difference of the QMC and LDA valence Compton 
profiles for lithium in the [100], [110] and [111] directions. We 
show the VMC and DMC results for the cell with 250 atoms 
and the VMC results for the cell with 686 atom. 



If we assume a constant valence density corresponding to 
r s = 3.25 a.u. and employ the QMC momentum distri- 
bution we computed for the electron gas, the resulting 
Lam-Platzman corrections are in qualitative agreement 
with the QMC correlation corrections of Fig |t|. Conse- 
quently, given the large difference between the experi- 
mental results by Sakurai et al. and the QMC profiles, 
the Lam-Platzman correction offer a satisfactory descrip- 
tion of electronic correlation in QMC even though they 
are isotropic and cannot resolve the directional differ- 
ences observed within QMC. 

In disagreement with the GW calculations Q , the QMC 
results indicate that electronic correlation in lithium only 
accounts for about 30% of the discrepancy between ex- 
perimental data and conventional band theoretical re- 
sults. A similar conclusion has been recently obtained in 
a comparison of QMC and experimental Compton pro- 
files for bulk silicon pig] . Other effects such as temper- 
ature ]ll[ ] or the interpretation of experimental results 
(in particular, the validity of the impulse approxima- 
tion 0|| ) may explain the difference between experi- 
ments and QMC calculations. Our calculations also show 
that QMC can not only offer a qualitative description of 



correlation effects in Compton profiles but, if finite size 
errors are properly taken care of, also resolve directional 
differences in the correlation corrections to the LDA re- 
sults and provide useful comparisons for new Compton 
scattering experiments Q . 
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APPENDIX A: 



LINEAR TETRAHEDRON 
METHOD 



To determine the momentum distribution within Kohn- 
Sham density functional theory, the Kohn-Sham valence 
wave functions and eigenvalues are evaluated on a very 
fine grid of ^-points within the irreducible wedge of the 
Brillouin zone. The Fourier coefficients of the orbitals 
with the corresponding eigenvalues are unfolded from the 
irreducible zone to full space and the momentum distri- 
bution is obtained by simply squaring the Fourier compo- 
nents. This procedure yields a square mesh with the mo- 
mentum distribution and single-particle energy defined 
at each point. 

To obtain the Compton profiles in a given direction, 
we integrate the momentum distribution using the lin- 
ear tetrahedron method as described by Lehmann and 
Taut J25| . The expression for the Compton profile (Eq. ||) 
is rewritten as 



J(P)= [ dpn(p)5(q-p -p)9(e(p) - e P ), 



(Al) 



where n(p) and e(p) are the Kohn-Sham single-particle 
momentum distribution and energy and e F the Fermi en- 
ergy. The ^-function is introduced to ensure that only 
states below the Fermi energy are included. If n(p) and 
e(p) have been computed on a square grid in momentum 
space, the mesh naturally divides space in cubes whose 
corners are defined by the grid points. Each cube is then 
divided in six tetrahedra and each tetrahedron is consid- 
ered in turn. Within each tetrahedron, n(p) and e(p) 
are linearly interpolated. The step function 0(e(p) — e F ) 
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may restrict the integration to only part of the tetrahe- 
dron: the primary tetrahedron is divided in secondary 
ones delimited by the boundaries e(p) = e F . For a given 
value of p, the surface of constant projection, q- p = p, is 
determined within each of the secondary tetrahedra and 
the contribution to J(p) is computed as the integral over 
this surface of the linear interpolation of the momentum 
distribution. 

Before computing the Compton profiles, the Fermi en- 
ergy e F must be estimated. For a given value of the Fermi 
energy, we apply the linear tetrahedron construction to 
determine the volume delimited by the Fermi surface and 
iteratively change e F so that 

Jdp6(e(p)-e F ) = ^- P l (A2) 



and increasing the plane-wave cutoff to 400 Ry has a 
negligible effect on the valence Compton profiles. 

TABLE I. Total energy of full-core bcc lithium versus the 
energy cutoff employed in the plane-wave calculation. The 
energy of full-core atomic lithium is -Eatom = —14.6682 Ry. 
All energies are in Rydberg. 



E cu t s 


-Etotal 


Els 


£2s 


300 


-14.6277 


-2.7442 


-0.3102 


400 


-14.6886 


-2.7619 


-0.3086 


1500 


-14.7985 


-2.7948 


-0.3058 


1800 


-14.8032 


-2.7962 


-0.3057 



APPENDIX C: QUANTUM MONTE CARLO 

where p F is the Fermi momentum. METHODS 



APPENDIX B: FULL-CORE PLANE-WAVE LDA 
CALCULATIONS 



As mentioned in Section IV, due to the lack of correct or- 



thogonalization between core and valence, the pseudopo- 
tential momentum distribution is too low at high mo- 
menta and, consequently, too high at low momenta. To 
account for the correct oscillatory behavior of the valence 
wave functions, the valence orbitals are usually com- 
puted within traditional all-electron schemes (LAPW or 
LMTO) or, alternatively, the true valence wave functions 
are reconstructed from the pseudized orbitals |p6| , p7f . 
On the other hand, several calculations in the literature 
have shown the feasibility of plane-wave calculations us- 
ing the unscreened Coulomb potential for first-row el- 
ements [ ^9|j2^p9[ |. We therefore decided to follow this 
more straightforward route and adopt a plane-wave ba- 
sis also with the — 3/r potential of lithium. 

In table ||, we show the convergence in the total energy 
and the Is and 2s eigenvalues at the T point for full- 
core solid lithium as a function of plane- wave cutoff. For 
zone sampling to selfconsistency, we used 14 special k- 
points. To obtain an accuracy of a few mRy on either 
the total energy or the Is eigenvalue, plane- wave cutoffs 
higher than 1800 Ry are required. On the other hand, 
the convergence in valence properties such as the 2 s level 
is achieved at significantly lower cutoffs. 

Bellaiche and Kunc |28| obtained convergence in the 
structural properties of solid LiH with a plane- wave cutoff 
of the order of 200 Ry and the valence Compton profiles 
were in good agreement with the ones reconstructed by 
simply orthogonalizing the valence wave functions to a 
core atomic orbital. Similarly, we find that the valence 
orbitals of lithium are well converged at about 300 Ry 



The variance minimization method |2J],|21|] consists of the 
minimization of the variance of the local energy over a 
set of N c configurations {Ri} sampled from the square 
of the best wave function available before we start the 
optimization, 



'opt 



E 



mi(Ri) 



N 



v{Ri)/^2w(Ri). (CI) 



E 



gUCSS 



is a guess for the energy of the state we are inter- 
ested in and w(R, t ) = \^(Ri)/^ (Ri)\ 2 . We do not allow 
the ratio of the weights to the average weight to exceed 
a maximum value. 

We compute the expectation value of various operators 
both in variational and diffusion Monte Carlo. In VMC, 
configurations are sampled from \1/ 2 using Metropolis 
Monte Carlo method and the expectation value of a given 
operator O is obtained from 







VMC 



1 N 

-Y 



Q^{R l ) 



(C2) 



The transition matrix consists of a drift-diffusion step 
with a time-step optimized to minimize the auto- 
correlation time. In DMC, the imaginary-time evolution 
operator exp(— Ht) is used to project out the ground 
state from the trial wave function within the fixed-node 
and the short-time approximations ifjp]] . The time-step 
error coming from the short-time approximation is negli- 
gible but the fixed-node error limits the accuracy of the 
results we obtain. 
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